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Artemisinin constitutes the frontline treatment to aid rapid clearance of parasitaemia and quick resolution 
of malarial symptoms. However, the widespread promiscuity about its mechanism of action is baffling. 
There is no consensus about the biochemical target of artemisinin but recent studies implicate haem and 
PfATP6 (a calcium pump). We investigated the role of iron and artemisinin on PfATP6, in search of a 
plausible mechanism of action, via density functional theory calculations, docking and molecular dynamics 
simulations. Results suggest that artemisinin gets activated by iron which in turn inhibits PfATP6 by closing 
the phosphorylation, nucleotide binding and actuator domains leading to loss of function of PfATP6 of the 
parasite and its death. The mechanism elucidated here should help in the design of novel antimalarials. 

Plasmodium falciparum is among the most prevalent species infecting humans killing more than two million 
people every year 1 ' 2 . Several reports have appeared in the recent past highlighting the complexity of its 
genomic organization, the metabolic pathways involved in infection and the potential for cure 3 " 6 . Resistance 
of these parasites to the traditional treatments has led to extensive work in discovering newer artemisinin analogs 
and derivatives 7,8 . Artemisinin (Fig. 1) is a sesquiterpene lactone endoperoxide containing a structural feature 
called peroxide bridge which is believed to be the key to its mode of action 9 . In spite of increasing popularity in the 
use of artemisinin based therapies, the mechanism of action of these sesquiterpene lactone endoperoxides has 
eluded researchers due to its controversial nature 10 . 

According to one school of thought, the cleavage of the peroxide bridge in presence of ferrous ion (Fe 2+ ) from 
haem forms highly reactive free radicals which rapidly rearrange to more stable carbon- centered radicals 1112 . It 
has been suggested that these artemisinin- derived free radicals chemically modify and inhibit a variety of parasite 
molecules, resulting in parasite's death 13 . A rich source of intracellular Fe 2+ is haem an essential component of 
hemoglobin, the malarial parasite is rich in haem iron derived from a breakdown of the host cell hemoglobin. It 
has long been suspected that Fe 2+ -haem is responsible for activating artemisinin inside the parasite 14 . Note that 
endoperoxides are known to be unstable, especially in the presence of iron. Theoretical studies on the cleavage of 
the peroxide bond by Taranto et al. 15 and Araujo et al. 16 suggest a thermodynamically favorable interaction 
between artemisinin and haem. 

Alternative views propose a direct interaction between artemisinin and PfATP6, a calcium pump. PfATP6 is 
the only SERCA-type Ca 2+ -ATPase present in the malarial parasite 17 . Some studies suggested that mutations on 
PfATP6 of the malarial parasite were able to modulate the affinity of artemisinin for the protein 1819 . The malarial 
enzyme was introduced into frog eggs and its activity was assessed, in the presence of artemisinin as well as 
thapsigargin, a known inhibitor for mammalian SERCA. It was found that both artemisinin and thapsigargin 
indeed inhibited SERCA. Furthermore, thapsigargin was found to interfere with the action of artemisinin 
indicating that they operated on malarial parasite through a similar mechanism 18 . In the parasite, the endoplasmic 
reticulum is situated outside the food vacuole throughout the cytoplasm, the same location as the proposed site of 
artemisinin action 14 . The parasite ingests and degrades up to 80% of the host-cell hemoglobin during its growth 
and replication process, in a compartment called food vacuole. Inside the malarial parasite, artemisinin, is 
activated by free iron neighboring PfATP6 in the endoplasmic reticulum. Abundance of ferrous ions within this 
vacuole, catalyze the cleavage of the peroxide bridge forming a highly reactive free radical 14 . 
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Figure 1 | Molecular structural formula of artemisinin. 

Modeling of PfATP6 enzyme based on templates HWO 20 ' 21 and 
2DQS 22 and docking of artemisinin derivative on to it was reported 
earlier. Jung et al. 20 and Naik et al. 21 obtained good correlation 
between their computed docking scores and in vitro antimalarial 
activities. However, no such correlations for several other antimalar- 
ial compounds with PfATP6 were observed 22 . Both 1IWO and 2DQS 
templates are from the same organism Oryctolagus Cunniculus 
bound with thapsigargin but having closed conformation of calcium 
pump. Reason for this contradictory results could be docking of 
artemisinin derivative to a static model of the enzyme or the starting 
structure which is in a closed conformation. 

In a nutshell, two types of mechanism of action for artemisinin 
have been proposed, one involving iron (or haem) 12,13 and the other 
independent of iron but involving PfATP6 18 . There have been several 
supportive 19,23 ' 24 and dismissive reports for these mechanisms 25 ' 26 
which makes the present study challenging 27 . It seems that a singular 
role of hemoglobin or PfATP6 may be insufficient to kill the parasite 
by artemisinin. We propose a mechanism of action of artemisinin 
which involves both iron and the binding of artemisinin to PfATP6 
enzyme. Understanding the issues that govern PfATP6 domain 
mobility, studied by molecular dynamics simulation and docking 
in the present study, may have profound implications for elucidating 
the detailed mechanism of binding of drug to this enzyme and in the 
design of new therapeutic agents, such as allosteric inhibitors, inten- 
ded to interfere with topological changes in the protein domains and 
thereby enzyme function. 

Amino acid sequence alignment clearly shows that PfATP6 con- 
serves all of the motifs and residues that are important for the struc- 
ture or function of a mammalian SERCA pump 28 . Both PfATP6 and 
mammalian SERCA shares 51% identity. It is also reported that 
PfATP6 displays structure and function similar to mammalian 
SERCA present in skeletal muscle 29 . Krishna et al. 30 and Salas- 
Burgos et al. 31 modeled the structure of PfATP6, (PDB ID: 1U5N) 
for which 1SU4 was used as a template which is an open conforma- 
tion crystal structure of calcium pump from the same organism 
Oryctolagus Cuniculus as studied earlier 20 " 22 except that previous 
studies used a closed conformation. In order to get accurate binding 
affinities upon docking, the enzyme has to undergo significant con- 
formational changes which no docking procedure allows. Closed 
conformation of PfATP6 models were built using II WO as template. 
Binding affinity computed using several docking and scoring meth- 
ods could not differentiate between inhibitor and non-inhibitor 
(such as deoxy-artemisinin) in closed model of PfATP6. However, 
our goal is to seek a plausible mechanism action of artemisinin bind- 
ing the enzyme whether that leads to similar change in conformation 
of PfATP6 as thapsigargin binding to mammalian SERCA (Open to 
closed conformations transformation leading to calcium pumping). 
Hence we studied from the model which is built using open con- 
formation template (1SU4) of the enzyme to see the dynamic effect of 



enzyme ligand complex along with the solvent. Alignment and 
superposition of the structures of the two proteins (1SU4 mam- 
malian SERCA and 1U5N (model of PfATP6) is shown in supple- 
mentary Fig. SI. In summary, PfATP6 retains all the characteristic 
features of the P-type ATPase family of SERCA pumps and can occur 
in broadly two conformations closed or open. There are 41 reported 
mammalian SERCA crystal structures in the RCSB 32 till date, most 
recent ones being 3W5A 33 and 4H1W 34 . All the structures except 
1SU4 and 2C9M are bound with different ligands and have closed 
or semi closed (El and E2) 35 conformations. A superposition of all 
these structures is shown in supplementary Fig. S2. The calculated 
average RMSD between the open conformation structure 1SU4 and 
the 39 closed conformations is ~13A Individual RMSDs of all these 
above mentioned proteins with respect to 1SU4 are compiled in 
supplementary table T2. 

Structure of PfATP6 is divided into three major parts: the cyto- 
plasmic headpiece, the stalk domain and the transmembrane domain 
(supplementary Fig. S3). The headpiece is composed of mainly the 
phosphorylation (P) domain, the nucleotide binding (N) domain, the 
hinge domain, and the actuator (A) domain or C-terminal region. 
The ligand binding site of PfATP6 was identified with thapsigargin as 
a reference ligand. The binding site is near the center of four trans- 
membrane helices. Small changes in these helices could lead to dis- 
ruption of the Ca 2+ binding sites and to alterations in the movement 
of ions to the cytosol or to luminal spaces 36 " 38 . 

Thapsigargin is known to induce structural change i.e. convert the 
mammalian SERCA from and open conformation to closed form. It 
is thus likely that the previous studies on artemisinin which used 
closed conformation were indecisive with respect to the activity as 
well as mechanism of action. We thus hypothesize that artemisinin 
or its Fe-adduct may act in manner similar to thapsigargin on mam- 
malian SERCA. It is thus of interest to study the open conformation 
of the enzyme and the influence of artemisinin and Fe- artemisinin 
adduct on the open conformation. Hence we adopted extensive 
molecular dynamics simulations of enzyme to analyze all the possible 
dynamic structures and intra domain motions in PfATP6. 

In the present work, we have carried out density functional theory 
calculations on artemisinin, haem and their complexes and con- 
ducted molecular docking experiments of these compounds with 
PfATP6. We have further carried out all atom molecular dynamics 
simulations for 100 ns on (i) Ca 2+ bound PfATP6, (ii) artemisinin 
bound PfATP6, and (iii) Fe-artemisinin adduct bound PfATP6 here 
after referred to as serca, arte-serca and Fe-arte-serca respectively. As 
a further control, we have also carried out molecular dynamics simu- 
lations on (iv) thapsigargin and mammalian SERCA complex. 

Results 

Density functional theory calculations on the complex between the 
haem group and artemisinin provided interaction energies, elec- 
tronic states, partial charges on each atom and geometry of the com- 
plex. Resultant structure of the complex was used for docking with 
PfATP6. The docked conformations of artemisinin and Fe-artemi- 
sinin adduct in the active site of PfATP6 are depicted in supplement- 
ary Fig S4. The binding pocket of PfATP6 consists of LEU 263, PHE 
264, GLN 267, ILE 977, ILE 981, ALA 985, ASN 1039, LEU 1040, ILE 
1041 and ASN 1042 residues. Interaction of artemisinin with the 
residue in the active site in the docked structure are consistent with 
the results reported earlier 20 22 . 

The essential massage of molecular dynamics (MD) simulations 
are captured in Fig. 2 in which systems (i) and (ii) viz. PfATP6 and 
artemisinin-PfATP6 complex show no major conformational 
changes. The two domains N & A, remain (—25 A) away from each 
other as in the native structure. However, they exhibit significant 
local fluctuations and some rearrangements. System (iii) namely, 
Fe-artemisinin adduct bound to PfATP6 shows a dramatic large 
conformational change resulting in a closure of N- and A- domains 
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Figure 2 | MD snapshots generated from individual trajectories after 
every 25 ns for PfATP6 enzyme (top row), artemisinin bound PfATP6 
(second row) and Fe-artemisinin adduct bound PfATP6 (third row). 

Nucleotide binding (N) domain and the actuator (A) domains are seen to 
be closing in the case of Fe-artemisinin adduct bound PfATP6 system. 



(see Fig. 2 and supplementary Fig. S5). Systerm (iv) too, namely 
thapsigargin bound to mammalian SERCA shows a similar large 
conformational change consistent with the experiment 37 (also see 
supplementary Fig. S6 and S7). 

Binding energies of the enzyme and ligand complex were calcu- 
lated 33 as average over 50 structures extracted after every two ns from 
the 100 ns simulation trajectories for the artemisinin-PfATP6, Fe- 
artemisinin-PfATP6 and thapsigargin-PfATP6 complexes and 
reported in Table 1. Fe- Artemisinin binds relatively strongly to 
PfATP6 but has lesser affinity towards mammalian SERCA 
(Table 1). However thapsigargin binds strongly with mammalian 
SERCA and shows reasonable affinity towards PfATP6. 
Experiments have also shown that thapsigargin interfered with the 
action of artemisinin 18 . It is clear from Table 1 that artemisinin and 
Fe activated artemisinin bind preferentially to PfATP6. thapsigargin 
which binds strongly in the active site of mammalian SERCA 36 has a 
larger volume 598 A 3 compared to artemisinin (258 A 3 ) as the active 
site of mammalian SERCA is large to fit thapsigargin tightly. 

To compare the flexibility of the structures, root mean square 
fluctuations were calculated (supplementary Fig. S8). Unbound pro- 
tein (SERCA) flexibility (supplementary Fig. S8a) was mostly 



localized in the loop regions of cytoplasmic domains. Differential 
fluctuation map of artemisinin-PfATP6 (Fig. S8b) Fe-artemisinin 
bound system (supplementary Fig. S8c) indicates suppressed and 
aggravated flexibility in specific regions. Fluctuations in the actuator 
domain, barring the residues which are connected to Ml helix, most 
of the residues display comparable fluctuations with respect to the 
native PfATP6 system. The nucleotide binding domain residues 
from 700 to 750 exhibit large fluctuations (supplementary Fig. S9). 
A comparison of arte-serca and Fe-arte-serca reveals the require- 
ment for activation of artemisinin by Fe as essential to bringing 
about the large conformational changes in PfATP6 leading to its 
inactivation. 

It is to be noted that (Fig. S3) the ATP-binding site in the N- 
domain is more than 25A away from the critical phosphorylation 
site SER 357, and about 80 A away from the Ca 2+ -binding site in both 
PfATP6 and mammalian SERCA. On the other hand, binding of 
Ca 2+ opens up the cleft in those domains. Thapsigargin, an inhibitor 
of mammalian SERCA, closes the N- and A- domain cleft, thus 
affecting the functioning of the enzyme 37 ' 39 . Interestingly, mam- 
malian SERCA, when inhibited by thapsigargin, leads to substantial 
conformational changes in the cytoplasmic headpiece. Three cyto- 
plasmic domains that were widely separated in the open structure of 
the Ca 2+ -bound structure were found to come close to each other in 
the Ca 2+ -free enzyme upon inhibition by thapsigargin 29,38 ' 39 . A sim- 
ilar conformational change is observed with Fe-artemisinin PfATP6 
complex but not with artemisinin PfATP6 or unbound PfATP6, 
Analogous spontaneous opening and reclosing of flaps of HIV-1 
protease with molecular dynamics simulations has been reported 
earlier 40 . 

An analysis of the MD trajectories shows that Fe-arte-serca during 
simulation moves inside the transmembrane region (supplementary 
Figs. S10 and SI 1) forming hydrogen bonds with ASN 1040, ILE 1039 
which are a part of the M5 helix. Additionally, bending of M5 helix is 
seen to occur accompanying a large conformational change in the 
open structure formed by the A-, P- and N-domains (Fig. 2). In Fe- 
arte-serca, changes in the topology of the enzyme structure are 
accompanied by a new interacting residue environment (Fig. 3a 
and 3b). Fe-artemisinin adduct moves towards the helical region 
forming hydrogen bond interactions with ASN 1039 and GLN 267 
allowing the ring part of artemisinin to interact favorably with LEU 
263, ILE 1041, ILE 1046, ILE 977 (Fig. 3a and 3b). 

The cleft formed by the P- and N-domains closes which could 
facilitate phosphorylation of SER 357 37 . Inter- atomic distances of 
residues from the nucleotide binding domain and the actuator 
domain decrease during the simulation of Fe-arte-serca (Fig. 4a). 
To obtain a quantitative view, we plotted the areas (Fig. 4b) of A-, 
H- and N-domains which is indicative of the movement of cytoplas- 
mic domains towards each other. 

To further understand the structural information with ther- 
modynamics of the open to closed transition, we constructed free 
energy landscapes from the MD trajectories performing principal 
component analysis (PC A). PC A 41,42 uses the actual dynamics of 
the protein to generate the motions of the native state of the protein 
conformation. In Fig. 5 the energy landscape in Fe-arte-serca system 



Table 1 Average binding free energies computed from the MD trajectories 

Molecule Target/Enzyme 


Binding Free Energy (kcal/mol) 


Artemisinin 


PfATP6* 


-6.5 


Fe-artemisinin 


PfATP6* 


-8.3 


Thapsigargin 


PfATP6* 


-6.7 


Artemisinin 


Mammalian SERCA (1SU4) 


-4.2 


Fe-artemisinin 


Mammalian SERCA (1SU4) 


-5.1 


Thapsigargin 


Mammalian SERCA (1SU4) 


-9.1 


'Modelled using 1 SU4 as described in text. 
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LEU 1046 



LEU 1041 




Figure 3 | Residues of PfATP6 in the vicinity of 5 A after 40 ns from simulation trajectories of (a) artemisinin (system (ii)) and (b) Fe-artemisinin 
adduct (system (iii)). Images were generated using VMD 56 



reveals characteristics that are not observed either with systems serca 
or arte-serca. There is a well-defined minimum structurally repre- 
senting the closed conformation state in Fe-arte-serca. The red 
region in all the plots has the least Gibbs free energy and their cor- 
responding structures, which are most probable, are pointed with 
arrows. 

A closer look (supplementary Fig. S4) at the ligand binding region 
shows that, Fe-artemisinin adduct forming directional hydrogen 
bonds with the backbone atoms of the functionally important resi- 
dues LEU 263 19 , ILE 977 20 and ASN 1039 20 . These H-bonds may 
guide Fe-artemisinin adduct deep into the membrane disrupting 
the helical region leading to strong and stable electrostatic inter- 
action between iron and side chain of ASN 980 from —10 A to 
~3A (supplementary Fig. S12). ASN 1039 (equivalent to mammalian 
PRO 827) and LEU263 (equivalent to mammalian GLU 255) seem to 
be giving selectivity to the artemisinin binding to PfATP6. 

A key question is how the ligand binding site and the phosphor- 
ylation site, 25A apart, communicate with each other. In this regard, 
it is important to note that the cytoplasmic end of M5 is integrated 
into the P-domain near the phosphorylation site and hydrogen 



bonded to the M4 helix. We analyzed the hydrogen bond linkage 
from the simulation trajectory of Fe-arte-serca system. There is an 
apparent linkage of network between Fe-artemisinin adduct and SER 
357 (phosphorylation site) which plays a pivotal role in the function- 
ing of SERCA. Also, SER 357 is directly connected to the nucleotide 
domain through a loop. Fe-artemisinin adduct makes hydrogen 
bond with ASN 1039 which is connected to SER 357 (phosphoryla- 
tion site) through hydrogen bonds, loop region and sheets. (Fig. 6). A 
tilt and lateral shift of M5 which leads to the large movement of 
nucleotide domain can also be attributed to the analogous linkage 
of Fe-artemisinin adduct with ILE 977 mentioned earlier. Main chain 
and H-bond connectivity are shown in Fig. 6a. Residue ASN 951, at 
one terminal of the M5 helix shows strong interaction with THR 362 
of phosphorylation domain, thus directly modulating the motion of 
the nucleotide domain. In addition to this, interaction of Fe-artemi- 
sinin adduct and residues of M4 (ASN 1039 and loop) helix also 
strongly influence the movement of the nucleotide domain. A net- 
work of hydrogen bonds between Fe-artemisinin adduct, LEU 263 
and LYS 259, attached to the actuator domain through a loop, reg- 
ulates the motion of the actuator domain as shown in Fig. 6b. There 




25 50 75 

Time (nanoseconds j 



25 50 75 

Time (nanoseconds) 



100 



a b 

Figure 4 | (a) Inter-atomic distances of backbone atom of VAL 221 of actuator domain and PHE 535 of nucleotide domain for (i) PfATP6 (red), (ii) 
artemisinin bound PfATP6 (green) and (iii) Fe-artemisinin adduct bound PfATP6 (blue) depicting the two residues coming closer from 25A to almost 8 
A for most part of the simulation (iii). However smaller variation in the distances were noticed with P/ATP6 (red), (26 to 32 A). For artemisinin bound 
PfATP6 system (green) (26- 20 A), (b) Area between the geometric centers of N-, H- and A- domains of PfATP6 system (red), artemisinin-PfATP6 system 
(green) and Fe-artemisinin adduct PfATP6 system (blue). Areas are calculated from triangles formed by the geometric centers of the three corresponding 
domains. 
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are other networks of ASP 358 alongside M5 and M4 where Fe- 
artemisinin adduct is located. However, the interactions mentioned 
above were found to be stronger indicated by a high survival of 
hydrogen bonds throughout the trajectory (almost 70% of the simu- 
lation run time taking hydrogen bond distance cut off as 3.2 A). These 
hydrogen bond linkages were not found with artemisinin bound 
complex during the simulation (supplementary Fig. SI 3). It is quite 
probable that many H-bond networks amongst P domain and the M5 
and M4 helices in association with other sites (Fig. 6) are responsible 
for the allosteric effect between the ligand binding site and the cyto- 
solic regions, causing the closure of two domains leading to inaccess- 
ibility of calcium ion at the calcium binding site. Thus, the events that 
occur at either site can be mechanically transmitted to the other site. 



Discussion 

Present study shows that Fe-artemisinin adduct upon binding to an 
open jaw like structure of PfATP6 comprising phosphorylation, nuc- 
leotide binding and actuator domains, is seen to cause large confor- 
mational changes in PfATP6 leading to a closure of these three 
domains causing inaccessibility of the calcium binding site and hence 
loss of function of Ca 2+ pump of the parasite and its resultant death. 
The dynamic behavior in terms of flexibility, conformation and the 
inhibitor- enzyme interaction of artemisinin and Fe-artemisinin 
adduct in the PfATP6 binding pocket are successfully explained by 
molecular simulations. The studies lead to the inference that Fe- 
artemisinin adduct inhibits PfATP6 through an allosteric mech- 
anism. Present studies on the binding mechanism of PfATP6 and 
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Fe-artemisinin adduct could prove helpful in the design of new anti- 
malarials. 

Methods 

Density functional theory (DFT) calculations. The geometries of both artemisinin 
and the haem group were fully optimized using density- functional theory as 
implemented in the ADF package 43 and as carried out by Araujo et al. The states with 
singlet, triplet and quintet multiplicities (S = 0, 1 and 2, respectively) were calculated 
for haem, while triplet, quintet, and septet states (S = 1,2, and 3, respectively) were 
calculated for the complex between haem and artemisinin. Triple zeta with single 
polarization (TZP) basis set and PBE functional 44 was used. The computed charges on 
each oxygen atoms are provided in Supplementary Table Tl. It was found that for 
interaction at short distances, due to artemisinin reduction by the haem group, the 
most stable complex has a septet spin state. Further computational details and the 
more insight into the interaction of haem and artemisinin can be found in that work. 
The resulting geometry which were of the same order as in the work of Araujo et al. 
were used for docking. 

Docking. Docking of artemisinin and Fe-artemisinin adduct to PfATP6 was carried 
out using Pardock 45 " 49 . Scoring of these docked structures was done by Bappl 49 . 
Details of the docking methodology are provided in the supplementary information. 
These docked structures were taken as input for molecular dynamics simulations. 

Molecular dynamics. All MD simulations were carried out using GROMACS 50 . 
Amber FF99SB 51 force field was used. The complexes obtained from Pardock were 
solvated in a box of TIP3P 52 water molecules (and 30 counter ions to ensure electro- 
neutrality) and minimized for an additional 2500 steps (1000SD + 1500CG). All 
these systems sizes comprised around —275,000 atoms. MD Simulations were started 
by slowly heating the solvent to 300 K over a period of 100 ps keeping the complex 
fixed. After this, the entire system was gradually brought to 300 K over a period of 
300 ps. The simulation was then carried out under NPT conditions for 100 ns. A 2 
femto second time step was used for integrating the equations of motion. Periodic 
boundary conditions were applied throughout the MD simulations, along with PME 
summation 53 for treating electrostatics. The temperature was kept constant by 
coupling heat bath through the Berendsen algorithm 54 using separate solute and 
solvent scaling. Pressure was adjusted by isotropic position scaling using a Berendsen 
like algorithm. Covalent bonds to hydrogen atoms were constrained by the SHAKE 
algorithm 55 . Convergence of energy, density were monitored. The initial 10 ns was 
treated as equilibration and the subsequent 100 ns was treated as production phase 
for structural and energy analyses. 
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